clear
set more off

* All required add-ons are stored in the /packages and /auxiliary folders
adopath ++ "$Wellness_WhatDoesWWDo/scripts/packages"
adopath ++ "$Wellness_WhatDoesWWDo/scripts/auxiliary"
mata: mata mlib index

*****************************
*** Set Color/Font Scheme ***
*****************************

local dkorange = `" "232 74 39" "'
local ltorange = `" "244 165 147" "'
local dkblue = `" "19 41 75" "'
local ltblue = `" "137 148 165" "'
local dkgray = `" "gs12" "'
local ltgray = `" "gs12*.5" "'
local vltgray = `" "gs15" "'

graph set window fontface "Times New Roman"

***********************
*** Create Figure 1 ***
***********************

* Figure 1 is created separately by "5b. make_figure1.do"
* This do file is called in "_run_all.do"
* The do file creates experimental_design.tex in /results/figures
* LaTeX needs to be run on that .tex file to create a PDF of the figure

*********************************************************
*** Simple Participation Rates among Treatment Groups ***
*********************************************************

use "$Wellness_WhatDoesWWDo/data/proc/wellness_analysis.dta", clear

gen participation = .
gen group = .
gen x_par = .
gen par_lab = ""

local j = 0

foreach v in screening_c hra_c activity_f_c activity_s_c screening2017_c hra_c_yr2 activity_f_c_yr2 activity_s_c_yr2 {
	
	summ `v' if treat==1
	local ++j
	replace participation = r(mean) in `j'
	
}

replace x_par = 2 + (_n - 1)*3 if _n <= 8

replace par_lab = string(100*participation,"%2.1f") + "%" if _n <= 8

replace group = 1 if inlist(_n,1,2,5,6) /*screening*/
replace group = 2 if inlist(_n,3,7) /*fall*/
replace group = 3 if inlist(_n,4,8) /*spring*/

twoway (bar participation x_par if group == 1, color(`dkorange') fintensity(inten100) base(0)) ///
	(bar participation x_par if group == 2, color(`dkblue') fintensity(inten100) base(0)) ///
	(bar participation x_par if group == 3, color(`dkgray') fintensity(inten100) base(0)) ///
	(scatter participation x_par, msymbol(i) mlabcolor(black) mlabel(par_lab) mlabposition(12)), ///
	ylabel(0 "0%" .2 "20%" .4 "40%" .6 "60%", gmin gmax glcolor(`vltgray') angle(0) labsize(small) noticks) ///
	xlabel(2 `" " " "Screening" "' 5 `" " " "HRA" "' ///
			6.5 `" " " " " " " "{bf:Year 1}" "' ///
			8 `" "Fall" "Activity" "' 11 `" "Spring" "Activity" "' ///
			14 `" " " "Screening" "' 17 `" " " "HRA" "' ///
			18.5 `" " " " " " " "{bf:Year 2}" "' ///
			20 `" "Fall" "Activity" "' 23 `" "Spring" "Activity" "', labsize(small) noticks) /// 
	yscale(range(0 .65) noline)  xscale(range(1 23) noline) ysize(4.25) xsize(5.5) ///
	graphregion(color(white) margin(r+5)) legend(off) xtitle("") ytitle("") ///
	name(par, replace)
	
gr_edit .AddLine added_lines editor 8 4.8 62 4.8
gr_edit .AddLine added_lines editor 8 4.8 8 5.5
gr_edit .AddLine added_lines editor 62 4.8 62 5.5
gr_edit .AddLine added_lines editor 63.8 4.8 117.8 4.8
gr_edit .AddLine added_lines editor 63.8 4.8 63.8 5.5
gr_edit .AddLine added_lines editor 117.8 4.8 117.8 5.5
gr_edit .added_lines[1].style.editstyle linestyle(width(medthin) color(black)) editcopy
gr_edit .added_lines[2].style.editstyle linestyle(width(medthin) color(black)) editcopy
gr_edit .added_lines[3].style.editstyle linestyle(width(medthin) color(black)) editcopy
gr_edit .added_lines[4].style.editstyle linestyle(width(medthin) color(black)) editcopy
gr_edit .added_lines[5].style.editstyle linestyle(width(medthin) color(black)) editcopy
gr_edit .added_lines[6].style.editstyle linestyle(width(medthin) color(black)) editcopy
 
graph export "$Wellness_WhatDoesWWDo/results/figures/participation_all.pdf", replace

********************************************************
*** Difference in mean spending and nonzero spending ***
********************************************************

use "$Wellness_WhatDoesWWDo/data/proc/wellness_analysis.dta", clear
drop if mi(treat)

gen y = .
gen x = 2 + (_n-1)*2 if _n <= 2
gen y_lab = ""
gen ci_hi = .
gen ci_lo = .

foreach y in spend_0816_0717 nonzero_spend_0816_0717 {

	if "`y'" == "spend_0816_0717" { 
		local weights "[aw=covg_0816_0717]"
		local y_lab `" "$" + string(y,"%9.0fc") "'
		local y_ax `" 0 "$0" 250 "$250" 500 "$500" 750 "$750" "'
		local fname "spend"
	}
	else if "`y'" == "nonzero_spend_0816_0717" { 
		local weights ""
		local y_lab `" string(100*y,"%2.1f") + "%" "'
		local y_ax `" 0 "0%" .25 "25%" .5 "50%" .75 "75%" 1 "100%" "'
		local fname "nonzero"
	}
	
	reg `y' i.treat `weights', robust
	
	replace y = _b[_cons] in 1
	replace y = _b[_cons] + _b[1.treat] in 2
	replace ci_hi = y + invttail(e(df_r),.025)*_se[1.treat] in 2
	replace ci_lo = y - invttail(e(df_r),.025)*_se[1.treat] in 2
	
	count if !treat & e(sample)
	local x_ax = `"2 "Control (N=`=string(r(N),"%9.0gc")')" "'
	count if treat & e(sample)
	local x_ax = `" `x_ax' 4 "Treatment (N=`=string(r(N),"%9.0gc")')" "'
	
	replace y_lab = `y_lab'
	
	twoway (bar y x if _n == 1, color(`dkblue') fintensity(inten100) base(0)) ///
		(bar y x if _n == 2, color(`dkorange') fintensity(inten100) base(0)) ///
		(rcap ci_hi ci_lo x if _n == 2, lwidth(vthin) color(black)) ///
		(scatter y x if _n == 1, msymbol(i) mlabcolor(black) mlabel(y_lab) mlabposition(12) mlabsize(medlarge)) ///
		(scatter ci_hi x if _n == 2, msymbol(i) mlabcolor(black) mlabel(y_lab) mlabposition(12) mlabsize(medlarge)), ///
		ylabel(`y_ax', angle(0) grid gmin gmax glcolor(`vltgray') noticks labsize(medlarge)) ///
		xlabel(`x_ax', noticks labsize(medlarge)) ///
		xtitle("") ytitle("") legend(off) graphregion(color(white)) ///
		xscale(noline range(1 5)) yscale(noline) ///
		name(`y'_dif, replace)
	
	graph export "$Wellness_WhatDoesWWDo/results/figures/ITT_`fname'.pdf", replace

}

****************************************************************************
*** Histograms for baseline spending, baseline income, and post spending ***
****************************************************************************

use "$Wellness_WhatDoesWWDo/data/proc/wellness_analysis.dta", clear

gen y = .
gen ci_hi = .
gen ci_lo = .

preserve

foreach v in spend_0715_0716 salary_0616 spend_0816_0717 {

	* Baseline vars: compare screened to unscreened, so limit sample to people in the treatment group
	if "`v'"=="spend_0715_0716" {		
		keep if treat==1
		local x_title = "Monthly Spending in Pre-Period (Dollars)"
		local group "hra_c"
		local scale ""
		local y_lab `" 0 "0%" .05 "5%" .1 "10%" .15 "15%" .2 "20%" .25 "25%" .3 "30%" .35 "35%" "'
		local control "Not Screened"
		local treat "Screened"
		local text1 ".28 18"
		local text2 ".27 18"
		local fname "spend_pre"
	}

	else if "`v'"=="salary_0616" {	
		keep if treat==1
		local x_title = "Salary in Pre-Period (Thousands of Dollars)"
		local group "hra_c"
		local scale "/1000"
		local y_lab `" 0 "0%" .05 "5%" .1 "10%" .15 "15%" .2 "20%" .25 "25%" .3 "30%" .35 "35%" "'
		local control "Not Screened"
		local treat "Screened"
		local text1 ".28 15.3"
		local text2 ".27 15.3"
		local fname "salary_pre"
	}	

	* Post vars: compare teratment to control, so limit sample to people in the study
	else if "`v'"=="spend_0816_0717" {
		drop if StudyArm=="Not Randomized":StudyArm 
		local x_title = "Monthly Spending in Post-Period (Dollars)"		
		local group "treat"
		local scale ""
		local y_lab `" 0 "0%" .05 "5%" .1 "10%" .15 "15%" .2 "20%" .25 "25%" .3 "30%" "'
		local control "Control         "
		local treat "Treated   "
		local text1 ".24 18"
		local text2 ".232 18"
		local fname "spend_post"
	}
	else error 1
	
	drop if mi(`v')
	gen group = `group'==1

	* Kolmogorov-Smirnov test
	ksmirnov `v', by(group)
	local KS_p = r(p)
	local KS_p = "= " + string(r(p),"%9.3f")
	if r(p) < 0.001 local KS_p = "< 0.001"
	
	* Create bins for the histogram, according to quantiles of !group (i.e., quantiles of either non-participants or control group)
	qui summ `v' if !group, d
	local p25 = round(`r(p25)', .01)
	local p50 = round(`r(p50)', .01)
	local p75 = round(`r(p75)', .01)
	local p90 = round(`r(p90)', .01)
	local p95 = round(`r(p95)', .01)
	gen `v'_cut = 1 + (`v'>0) + (`v'>`r(p25)') + (`v'>`r(p50)') + (`v'>`r(p75)') + (`v'>`r(p90)') + (`v'>`r(p95)')
	
	if "`v'"=="salary_0616" {
		local x_ax `" 1.5 `"($0–`=round(`p25'`scale')']"' 6.5 "($`=round(`p25'`scale')'–`=round(`p50'`scale')']" 11.5 "($`=round(`p50'`scale')'–`=round(`p75'`scale')']" 16.5 "($`=round(`p75'`scale')'–`=round(`p90'`scale')']" 21.5 "($`=round(`p90'`scale')'–`=round(`p95'`scale')']" 26.5 "$`=round(`p95'`scale')'+" "'
	}
	else {
		local x_ax `" 1.5 "$0" 6.5 `"(0–`=round(`p25'`scale')']"' 11.5 "(`=round(`p25'`scale')'–`=round(`p50'`scale')']" 16.5 "(`=round(`p50'`scale')'–`=round(`p75'`scale')']" 21.5 "(`=round(`p75'`scale')'–`=round(`p90'`scale')']" 26.5 "(`=round(`p90'`scale')'–`=round(`p95'`scale')']" 31.5 "`=round(`p95'`scale')'+" "'
	}

	* calculate chi-squared test
	tab `v'_cut group, col nofreq chi2
	local cs_p = "= " + string(r(p),"%9.3f")
	if r(p) < 0.001 local cs_p = "< 0.001"
	
	* Bins and confidence intervals
	qui forvalues i = 1/7 {
	
		qui gen `v'_cut_`i' = `v'_cut==`i'

		qui _regress `v'_cut_`i' group, robust
		replace y = _b[_cons] if _n == 2*`i' - 1
		replace y = _b[_cons] + _b[group] if _n == 2*`i'
		replace ci_hi = y + invttail(e(df_r),.025)*_se[group] if _n == 2*`i'
		replace ci_lo = y - invttail(e(df_r),.025)*_se[group] if _n == 2*`i'
		
		drop `v'_cut_`i'
		
	}
	
	gen x = .5*(5*_n - 6 + mod(_n,2)*3) if _n <= 14
	
	if "`v'"=="salary_0616" {
	
		replace x = x - 5
		replace x = . if x < 0
	
	}
	
	twoway (bar y x if mod(_n,2), color(`dkblue') fintensity(inten100) base(0)) ///
		(bar y x if !mod(_n,2), color(`dkorange') fintensity(inten100) base(0)) ///
		(rcap ci_hi ci_lo x if !mod(_n,2), lwidth(vthin) color(black)), ///
		ylabel(`y_lab', angle(0) grid gmin gmax glcolor(`vltgray') noticks labsize(medsmall)) ///
		xlabel(`x_ax', noticks labsize(small)) ysize(4.25) xsize(5.5) ///
		xtitle("`x_title'", size(medium)) ytitle("") graphregion(color(white)) xscale(noline) yscale(noline) ///
		legend(order(1 "`control'" 2 "`treat'") symxsize(*.2) bmargin(medium) ring(0) bplacement(ne) size(small) region(color(white)) position(2)) ///
		text(`text1' "Pearson's chi-squared test for equality: p-value `cs_p'", place(e) size(vsmall)) ///
		text(`text2' "Kolmogorov-Smirnov test for equality: p-value `KS_p'", place(e) size(vsmall)) ///
		name(`v'_distribution, replace)
		
	graph export "$Wellness_WhatDoesWWDo/results/figures/histogram_`fname'.pdf", replace
	
	restore, preserve
}

**************************
*** Lit Review Figures ***
**************************

local dkorange = `" "232 74 39" "'
local ltorange = `" "244 165 147" "'
local dkblue = `" "19 41 75" "'
local ltblue = `" "137 148 165" "'
local dkgray = `" "gs12" "'
local ltgray = `" "gs12*.5" "'
local vltgray = `" "gs15" "'

graph set window fontface "Times New Roman"

local lasso_dir = "$Wellness_WhatDoesWWDo/data/proc/0816_0717/lasso"
local winsorize_dir = "$Wellness_WhatDoesWWDo/results/intermediate_files/0816_0717/winsorization/estimates"
local estimate_dir = "$Wellness_WhatDoesWWDo/results/intermediate_files/0816_0717/treat_effects/estimates"

* Make Data Set using excel sheet with details about each study

import excel using "$Wellness_WhatDoesWWDo/data/raw/lit_review/wellness_lit_review.xlsx", sheet("Export") firstrow clear

varabbrev {

rename Title title
rename Drop_Dup drop_dup
rename Review review
rename Source source
rename Data study_year
rename Outcome outcome
rename BaselineLevel baseline_level
rename OriginalEffect effect_level
rename percent_change effect_percent
rename Checked checked
rename Pre_Post pre_post
rename itt_tot estimate_type
rename itt_participation itt_share
rename OriginalTimeIntervalinmonth timelength
rename standarderror standard_error

}

keep title review timelength source study_year outcome baseline_level /// 
	effect_level effect_percent checked pre_post estimate_type study_type /// 
	estimate_type itt_share drop_dup alt_effect standard_error p_value ///
	z_score p_upper_bound has_se

drop if effect_percent == .

replace outcome = "health" if outcome == "health care costs" | outcome == "health utilization"
replace outcome = "absenteeism" if outcome == "absenteeism costs"
replace source = "Original Paper" if source != "Review Article"
replace drop_dup = 0 if drop_dup == .
replace drop_dup = 0 if mi(drop_dup)
destring timelength, replace

foreach x of varlist review source outcome study_type estimate_type {
	
	rename `x' temp_var
	encode temp_var, generate(`x')
	drop temp_var

}

** Create variables for publication bias analysis

gen pub_effect = effect_level
replace pub_effect = alt_effect if alt_effect != .

** Generate standard error variable, using se, p-value, or z-score

gen pub_se = standard_error

replace pub_se = abs(pub_effect/invnormal(1 - p_value/2)) ///
	if pub_se == . & p_value != .

replace pub_se = abs(pub_effect/z_score) if pub_se == . & p_value == . ///
	& z_score != .

order title review source study_year outcome baseline_level effect_level effect_percent checked pre_post study_type estimate_type itt_share drop_dup timelength pub_effect pub_se

sum checked, meanonly
assert r(mean) == 1 /* all estimates have been double-checked */

cap mkdir "$Wellness_WhatDoesWWDo/data/proc/lit_review"
save "$Wellness_WhatDoesWWDo/data/proc/lit_review/lit_review.dta", replace

* Graph Distribution of Estimates

** ITT Effects for health costs, 1% Winsorized

* Load estimates and calculate point estimate and CI

estimates use "`winsorize_dir'/ITT_Lasso_Winsorized_spendw01_0816_0717.ster"

use "`lasso_dir'/lasso_data_spend_0816_0717.dta", clear
* Merge in winsorized spending
merge 1:1 AnalysisID using "$Wellness_WhatDoesWWDo/data/proc/wellness_analysis.dta", assert(match using) keep(match) nogen keepusing(spendw*)
assert _N==4834

sum spendw01_0816_0717 [aweight = covg_0816_0717] if in_lasso_IV == 1 & treat == 0

local health_itt = _b[rhs]/r(mean)
local health_itt_lb = (_b[rhs] - invttail(e(df_r),0.025)*_se[rhs])/r(mean)
local health_itt_ub = (_b[rhs] + invttail(e(df_r),0.025)*_se[rhs])/r(mean)

* Load estimates from prior studies

use "$Wellness_WhatDoesWWDo/data/proc/lit_review/lit_review.dta", clear
drop if drop_dup
keep if estimate_type == "itt":estimate_type & outcome == "health":outcome

* Create an observation for our point estimate

set obs `=_N+1'
replace title = "Paper Estimate" in `=_N'
replace effect_percent = `health_itt' in `=_N'
gen main = _n == _N

* Put all point estimates into one of 50 bins

sum effect_percent, detail

local max_x = max(ceil(100*r(max))/100,`health_itt_ub')
local delta = (`max_x' + 1)/50
gen bin = .

forval x = 1/50 {

	replace bin = `max_x' - `x'*`delta' + .5*`delta' if effect_percent < `max_x' - `=`x'-1'*`delta' & effect_percent >= `max_x' - `x'*`delta'

}

egen bin_code = group(bin), label

di "`health_itt_lb' `health_itt' `health_itt_ub'"

* Label point estimates that are out(side) of 95% CI
* Do it by bin, then move estimates near the boundary to the next bin if they are "out"

bysort bin_code (effect_percent): gen out = effect_percent[1] > `health_itt_ub' | effect_percent[_N] < `health_itt_lb'

replace bin_code = bin_code - 1 if !out & effect_percent < `health_itt_lb'
replace bin_code = bin_code + 1 if !out & effect_percent > `health_itt_ub'
replace out = 1 if out == 0 & (effect_percent > `health_itt_ub' | effect_percent < `health_itt_lb')

rename bin bin2
decode bin_code, gen(bin)
destring bin, replace

gen order = .
replace order = 4 if main
replace order = 3 if !main & study_type == "RCT":study_type
replace order = 2 if !main & study_type != "RCT":study_type & pre_post
replace order = 1 if !main & study_type != "RCT":study_type & !pre_post

bysort bin_code (order), sort: gen counter = _N + 1 - _n

sum bin if main, meanonly
local itt_x = r(mean)
sum counter if main, meanonly
local itt_y = r(mean)

twoway scatter counter bin if !out & order == 3, mcolor(`ltblue') msymbol(S) || ///
scatter counter bin if !out & order == 2, mcolor(`ltblue') msymbol(T) || ///
scatter counter bin if !out & order == 1, mcolor(`ltblue') msymbol(O) || /// 
scatter counter bin if out & order == 3, mcolor(`ltblue') msymbol(Sh) || ///
scatter counter bin if out & order == 2, mcolor(`ltblue') msymbol(Th) || ///
scatter counter bin if out & order == 1, mcolor(`ltblue') msymbol(Oh) || ///
scatter counter bin if order == 4, msymbol(none) ||, ///
text(`itt_y' `itt_x' "★", color(`dkorange')) ///
yscale(off r(1 20)) ylabel(none,nogrid) graphregion(color(white)) ///
xlabel(-.5 "-50%" -.25 "-25%" 0 "0%" .25 "25%") xtitle("Percent Change in Medical Spending (ITT)") ytitle("") ///
legend(label(1 "RCT") label(2 "Pre/Post") label(3 "Other") ///
label(4 "RCT") label(5 "Pre/Post") label(6 "Other") ///
label(7 "RCT") ///
order(- "IL Estimates:" 7 - "95% Confidence Interval" - ""  /// 
- "Prior Studies:" - "★" - "[ ]" - "" ///
- "Not Ruled Out:" 1 2 3 /// 
- "Ruled Out:" 4 5 6) rows(4)) ///
name(health_itt, replace)

sum counter
local top = r(max) + 1
local lb = `health_itt_lb' - .01

gr_edit .plotregion1.AddLine added_lines editor `lb' .5 `lb' `top'
gr_edit .plotregion1.AddLine added_lines editor `lb' .5 `=`lb'+ 0.02' .5
gr_edit .plotregion1.AddLine added_lines editor `lb' `top' `=`lb'+ 0.02' `top'
gr_edit .plotregion1.AddLine added_lines editor `health_itt_ub' .5 `health_itt_ub' `top'
gr_edit .plotregion1.AddLine added_lines editor `health_itt_ub' .5 `=`health_itt_ub'- 0.02' .5
gr_edit .plotregion1.AddLine added_lines editor `health_itt_ub' `top' `=`health_itt_ub'- 0.02' `top'
gr_edit .plotregion1.added_lines[1].style.editstyle  linestyle( width(thin) color(red))
gr_edit .legend.plotregion1.key[4].DragBy 0 -26
gr_edit .legend.plotregion1.key[7].DragBy 0 -26
gr_edit .legend.plotregion1.label[12].DragBy 0 -26
gr_edit .legend.plotregion1.label[16].DragBy 0 -26
gr_edit .legend.plotregion1.label[9].DragBy 0 2.8
gr_edit .legend.plotregion1.label[13].DragBy .0 2.8
gr_edit .legend.plotregion1.label[6].DragBy 5.188841396576663 -4.425774722311327
gr_edit .legend.plotregion1.label[7].DragBy 5.188841396576663 -4.425774722311327
gr_edit .legend.plotregion1.label[6].style.editstyle color(`dkorange') editcopy
gr_edit .legend.plotregion1.label[7].style.editstyle color(red) editcopy
gr_edit .legend.plotregion1.label[6].style.editstyle size(small) editcopy
gr_edit .plotregion1.textbox1.style.editstyle vertical(top) editcopy
gr_edit .plotregion1.textbox1.style.editstyle size(small) editcopy

graph export "$Wellness_WhatDoesWWDo/results/figures/lit_review_health_ITT.pdf", replace

twoway scatter counter bin if !out & order == 3, mcolor(`ltblue') msymbol(S) || ///
scatter counter bin if !out & order == 2, mcolor(`ltblue') msymbol(T) || /// 
scatter counter bin if !out & order == 1, mcolor(`ltblue') msymbol(O) || /// 
scatter counter bin if out & order == 3, mcolor(`ltblue') msymbol(Sh) || ///
scatter counter bin if out & order == 2, mcolor(`ltblue') msymbol(Th) || ///
scatter counter bin if out & order == 1, mcolor(`ltblue') msymbol(Oh) || ///
scatter counter bin if order == 4, msymbol(none) ||, ///
text(`itt_y' `itt_x' "★", color(`dkorange')) ///
yscale(off r(1 20)) ylabel(none,nogrid) graphregion(color(white)) legend(off) ///
xlabel(-.5 "-50%" -.25 "-25%" 0 "0%" .25 "25%") xtitle("Percent Change in Medical Spending (ITT)") ytitle("") ///
name(health_itt2, replace)

sum counter
local top = r(max) + 1
local lb = `health_itt_lb' - .01

gr_edit .plotregion1.AddLine added_lines editor `lb' .65 `lb' `top'
gr_edit .plotregion1.AddLine added_lines editor `lb' .65 `=`lb'+ 0.02' .65
gr_edit .plotregion1.AddLine added_lines editor `lb' `top' `=`lb'+ 0.02' `top'
gr_edit .plotregion1.AddLine added_lines editor `health_itt_ub' .65 `health_itt_ub' `top'
gr_edit .plotregion1.AddLine added_lines editor `health_itt_ub' .65 `=`health_itt_ub'- 0.02' .65
gr_edit .plotregion1.AddLine added_lines editor `health_itt_ub' `top' `=`health_itt_ub'- 0.02' `top'
gr_edit .plotregion1.added_lines[1].style.editstyle  linestyle( width(thin) color(red))
gr_edit .plotregion1.textbox1.style.editstyle vertical(top) editcopy
gr_edit .plotregion1.textbox1.style.editstyle size(small) editcopy

graph export "$Wellness_WhatDoesWWDo/results/figures/lit_review_health_ITT_no_legend.pdf", replace

** IV Effects for health costs, 1% Winsorized

* Load estimates and calculate point estimate and CI

estimates use "`winsorize_dir'/IV_Lasso_Winsorized_spendw01_0816_0717.ster"

use "`lasso_dir'/lasso_data_spend_0816_0717.dta", clear

* JR edit
*sum spend_0816_0717 [aweight= covg_0816_0717] if in_lasso_IV == 1 & hra_c_nomiss == 1

* Merge in winsorized spending
merge 1:1 AnalysisID using "$Wellness_WhatDoesWWDo/data/proc/wellness_analysis.dta", assert(match using) keep(match) nogen keepusing(spendw*)
assert _N==4834

sum spendw01_0816_0717 [aweight= covg_0816_0717] if in_lasso_IV == 1 & hra_c_nomiss == 1

local health_iv = _b[rhs]/(r(mean) - _b[rhs])
local health_iv_lb = (_b[rhs] -  invnormal(0.975)*_se[rhs])/(r(mean) - _b[rhs])
local health_iv_ub = (_b[rhs] +  invnormal(0.975)*_se[rhs])/(r(mean) - _b[rhs])

* Load OLS results

estimates use "`estimate_dir'/OLS_Lasso_spend_0816_0717.ster"

use "`lasso_dir'/lasso_data_spend_0816_0717.dta", clear
sum spend_0816_0717 [aweight= covg_0816_0717] if in_lasso_OLS == 1 & hra_c_nomiss == 0

local health_ols = _b[rhs]/(r(mean) - _b[rhs])

* Load estimates from prior studies

use "$Wellness_WhatDoesWWDo/data/proc/lit_review/lit_review.dta", clear
drop if drop_dup
keep if estimate_type == "tot":estimate_type & outcome == "health":outcome

* Create an observation for our point estimate

set obs `=_N+2'
replace title = "Paper Estimate IV" in `=_N - 1'
replace title = "Paper Estimate OLS" in `=_N'
replace effect_percent = `health_iv' in `=_N - 1'
replace effect_percent = `health_ols' in `=_N'
gen main_iv = _n == _N - 1
gen main_ols = _n == _N
gen main = _n >= _N - 1

* Put all point estimates into one of 50 bins

sum effect_percent, detail

local max_x = max(ceil(100*r(max))/100,`health_iv_ub')
local delta = (`max_x' + 1)/50
gen bin = .

forval x = 1/50 {

	replace bin = `max_x' - `x'*`delta' + .5*`delta' if effect_percent < `max_x' - `=`x'-1'*`delta' & effect_percent >= `max_x' - `x'*`delta'

}

egen bin_code = group(bin), label

di "`health_iv_lb' `health_iv' `health_iv_ub'"

* Label point estimates that are out(side) of 95% CI
* Do it by bin, then move estimates near the boundary to the next bin if they are "out"

bysort bin_code (effect_percent): gen out = effect_percent[1] > `health_iv_ub' | effect_percent[_N] < `health_iv_lb'

replace bin_code = bin_code - 1 if !out & effect_percent < `health_iv_lb'
replace bin_code = bin_code + 1 if !out & effect_percent > `health_iv_ub'
replace out = 1 if out == 0 & (effect_percent > `health_iv_ub' | effect_percent < `health_iv_lb')

rename bin bin2
decode bin_code, gen(bin)
destring bin, replace

gen order = .
replace order = 5 if main_iv
replace order = 4 if main_ols
replace order = 3 if !main & study_type == "RCT":study_type
replace order = 2 if !main & study_type != "RCT":study_type & pre_post
replace order = 1 if !main & study_type != "RCT":study_type & !pre_post

bysort bin_code (order), sort: gen counter = _N + 1 - _n

sum bin if main_iv, meanonly
local iv_x = r(mean)
sum counter if main_iv, meanonly
local iv_y = r(mean)

sum bin if main_ols, meanonly
local ols_x = r(mean)
sum counter if main_ols, meanonly
local ols_y = r(mean)

twoway scatter counter bin if !out & order == 3, mcolor(`ltblue') msymbol(S) || ///
scatter counter bin if !out & order == 2, mcolor(`ltblue') msymbol(T) || /// 
scatter counter bin if !out & order == 1, mcolor(`ltblue') msymbol(O) || /// 
scatter counter bin if out & order == 3, mcolor(`ltblue') msymbol(Sh) || ///
scatter counter bin if out & order == 2, mcolor(`ltblue') msymbol(Th) || ///
scatter counter bin if out & order == 1, mcolor(`ltblue') msymbol(Oh) || ///
scatter counter bin if order == 4, msymbol(none) || ///
scatter counter bin if order == 5, msymbol(none) ||, ///
text(`iv_y' `iv_x' "★", color(`dkorange')) text(`ols_y' `ols_x' "☆", color(`dkorange')) ///
yscale(off r(1 20)) ylabel(none,nogrid) graphregion(color(white)) ///
xlabel(-1 "-100%" -.5 "-50%" 0 "0%" .5 "50%" 1 "100%" 1.5 "150%") xtitle("Percent Change in Medical Spending (TOT)") ytitle("") ///
legend(label(1 "RCT") label(2 "Pre/Post") label(3 "Other") ///
label(4 "RCT") label(5 "Pre/Post") label(6 "Other") ///
label(7 "OLS") label(8 "RCT") ///
order(- "IL Estimates:" 8 - "95% Confidence Interval" 7 /// 
- "Prior Studies:" - "★" - "[ ]" - "☆" ///
- "Not Ruled Out:" 1 2 3 /// 
- "Ruled Out:" 4 5 6) rows(4)) ///
name(health_iv, replace)

sum counter
local top = r(max) + 1
local lb = `health_iv_lb' - 0.005

gr_edit .plotregion1.AddLine added_lines editor `lb' .5 `lb' `top'
gr_edit .plotregion1.AddLine added_lines editor `lb' .5 `=`lb'+ 0.02' .5
gr_edit .plotregion1.AddLine added_lines editor `lb' `top' `=`lb'+ 0.02' `top'
gr_edit .plotregion1.AddLine added_lines editor `health_iv_ub' .5 `health_iv_ub' `top'
gr_edit .plotregion1.AddLine added_lines editor `health_iv_ub' .5 `=`health_iv_ub'- 0.02' .5
gr_edit .plotregion1.AddLine added_lines editor `health_iv_ub' `top' `=`health_iv_ub'- 0.02' `top'
gr_edit .plotregion1.added_lines[1].style.editstyle  linestyle( width(thin) color(red))
gr_edit .legend.plotregion1.label[6].DragBy 5.188841396576663 -4.425774722311327
gr_edit .legend.plotregion1.label[7].DragBy 5.188841396576663 -4.425774722311327
gr_edit .legend.plotregion1.label[8].DragBy 5.188841396576663 -4.425774722311327
gr_edit .legend.plotregion1.key[5].DragBy 0 -26
gr_edit .legend.plotregion1.key[8].DragBy 0 -26
gr_edit .legend.plotregion1.label[12].DragBy 0 -26
gr_edit .legend.plotregion1.label[16].DragBy 0 -26
gr_edit .legend.plotregion1.label[9].DragBy 0 2.8
gr_edit .legend.plotregion1.label[13].DragBy .0 2.8
gr_edit .legend.plotregion1.label[7].style.editstyle color(red) editcopy
gr_edit .legend.plotregion1.label[6].style.editstyle color(`dkorange') editcopy
gr_edit .legend.plotregion1.label[8].style.editstyle color(`dkorange') editcopy
gr_edit .legend.plotregion1.label[6].style.editstyle size(small) editcopy
gr_edit .legend.plotregion1.label[8].style.editstyle size(small) editcopy
gr_edit .plotregion1.textbox1.style.editstyle vertical(top) editcopy
gr_edit .plotregion1.textbox1.style.editstyle size(small) editcopy
gr_edit .plotregion1.textbox2.style.editstyle vertical(top) editcopy
gr_edit .plotregion1.textbox2.style.editstyle size(small) editcopy

graph export "$Wellness_WhatDoesWWDo/results/figures/lit_review_health_TOT.pdf", replace

twoway scatter counter bin if !out & order == 3, mcolor(`ltblue') msymbol(S) || ///
scatter counter bin if !out & order == 2, mcolor(`ltblue') msymbol(T) || ///
scatter counter bin if !out & order == 1, mcolor(`ltblue') msymbol(O) || /// 
scatter counter bin if out & order == 3, mcolor(`ltblue') msymbol(Sh) || ///
scatter counter bin if out & order == 2, mcolor(`ltblue') msymbol(Th) || ///
scatter counter bin if out & order == 1, mcolor(`ltblue') msymbol(Oh) || ///
scatter counter bin if order == 4, msymbol(none) || ///
scatter counter bin if order == 5, msymbol(none) ||, ///
text(`iv_y' `iv_x' "★", color(`dkorange')) text(`ols_y' `ols_x' "☆", color(`dkorange')) ///
yscale(off r(1 20)) ylabel(none,nogrid) graphregion(color(white)) legend(off) ///
xlabel(-1 "-100%" -.5 "-50%" 0 "0%" .5 "50%" 1 "100%" 1.5 "150%") xtitle("Percent Change in Medical Spending (TOT)") ytitle("") ///
name(health_iv2, replace)

sum counter
local top = r(max) + 1
local lb = `health_iv_lb' - 0.005

gr_edit .plotregion1.AddLine added_lines editor `lb' .65 `lb' `top'
gr_edit .plotregion1.AddLine added_lines editor `lb' .65 `=`lb'+ 0.02' .65
gr_edit .plotregion1.AddLine added_lines editor `lb' `top' `=`lb'+ 0.02' `top'
gr_edit .plotregion1.AddLine added_lines editor `health_iv_ub' .65 `health_iv_ub' `top'
gr_edit .plotregion1.AddLine added_lines editor `health_iv_ub' .65 `=`health_iv_ub'- 0.02' .65
gr_edit .plotregion1.AddLine added_lines editor `health_iv_ub' `top' `=`health_iv_ub'- 0.02' `top'
gr_edit .plotregion1.added_lines[1].style.editstyle  linestyle( width(thin) color(red))
gr_edit .plotregion1.textbox1.style.editstyle vertical(top) editcopy
gr_edit .plotregion1.textbox1.style.editstyle size(small) editcopy
gr_edit .plotregion1.textbox2.style.editstyle vertical(top) editcopy
gr_edit .plotregion1.textbox2.style.editstyle size(small) editcopy

graph export "$Wellness_WhatDoesWWDo/results/figures/lit_review_health_TOT_no_legend.pdf", replace

** ITT effects for absenteeism

* Load estimates and calculate point estimate and CI

estimates use "`estimate_dir'/ITT_Lasso_sickleave_0816_0717.ster"

use "`lasso_dir'/lasso_data_sickleave_0816_0717.dta", clear
sum sickleave_0816_0717 [aweight= sickdays_eligible_0816_0717] if in_lasso_IV == 1 & treat == 0

local absent_itt = _b[rhs]/r(mean)
local absent_itt_lb = (_b[rhs] - invttail(e(df_r),0.025)*_se[rhs])/r(mean)
local absent_itt_ub = (_b[rhs] + invttail(e(df_r),0.025)*_se[rhs])/r(mean)

* Load estimates from prior studies

use "$Wellness_WhatDoesWWDo/data/proc/lit_review/lit_review.dta", clear
drop if drop_dup
keep if estimate_type == "itt":estimate_type & outcome == "absenteeism":outcome

* Create an observation for our point estimate

set obs `=_N+1'
replace title = "Paper Estimate" in `=_N'
replace effect_percent = `absent_itt' in `=_N'
gen main = _n == _N

* Put all point estimates into one of 50 bins

sum effect_percent, detail

local max_x = max(ceil(100*r(max))/100,`absent_itt_ub')
local delta = (`max_x' + 1)/50
gen bin = .

forval x = 1/50 {

	replace bin = `max_x' - `x'*`delta' + .5*`delta' if effect_percent < `max_x' - `=`x'-1'*`delta' & effect_percent >= `max_x' - `x'*`delta'

}

egen bin_code = group(bin), label

di "`absent_itt_lb' `absent_itt' `absent_itt_ub'"

* Label point estimates that are out(side) of 95% CI
* Do it by bin, then move estimates near the boundary to the next bin if they are "out"

bysort bin_code (effect_percent): gen out = effect_percent[1] > `absent_itt_ub' | effect_percent[_N] < `absent_itt_lb'

replace bin_code = bin_code - 1 if !out & effect_percent < `absent_itt_lb'
replace bin_code = bin_code + 1 if !out & effect_percent > `absent_itt_ub'
replace out = 1 if out == 0 & (effect_percent > `absent_itt_ub' | effect_percent < `absent_itt_lb')

rename bin bin2
decode bin_code, gen(bin)
destring bin, replace

gen order = .
replace order = 4 if main
replace order = 3 if !main & study_type == "RCT":study_type
replace order = 2 if !main & study_type != "RCT":study_type & pre_post
replace order = 1 if !main & study_type != "RCT":study_type & !pre_post

bysort bin_code (order), sort: gen counter = _N + 1 - _n

sum bin if main, meanonly
local itt_x = r(mean)
sum counter if main, meanonly
local itt_y = r(mean)

twoway scatter counter bin if !out & order == 3, mcolor(`ltblue') msymbol(S) || ///
scatter counter bin if !out & order == 2, mcolor(`ltblue') msymbol(T) || /// 
scatter counter bin if !out & order == 1, mcolor(`ltblue') msymbol(O) || /// 
scatter counter bin if out & order == 3, mcolor(`ltblue') msymbol(Sh) || ///
scatter counter bin if out & order == 2, mcolor(`ltblue') msymbol(Th) || ///
scatter counter bin if out & order == 1, mcolor(`ltblue') msymbol(Oh) || ///
scatter counter bin if order == 4, msymbol(none) ||, ///
text(`itt_y' `itt_x' "★", color(`dkorange')) ///
yscale(off r(1 20)) ylabel(none,nogrid) graphregion(color(white)) ///
xlabel(-.75 "-75%" -.5 "-50%" -.25 "-25%" 0 "0%" .25 "25%") xtitle("Percent Change in Absenteeism (ITT)") ytitle("") ///
legend(label(1 "RCT") label(2 "Pre/Post") label(3 "Other") ///
label(4 "RCT") label(5 "Pre/Post") label(6 "Other") ///
label(7 "RCT") ///
order(- "IL Estimates:" 7 - "95% Confidence Interval" - ""  /// 
- "Prior Studies:" - "★" - "[ ]" - "" ///
- "Not Ruled Out:" 1 2 3 /// 
- "Ruled Out:" 4 5 6) rows(4)) ///
name(absent_itt, replace)

sum counter
local top = r(max) + 1
local lb = `absent_itt_lb'

gr_edit .plotregion1.AddLine added_lines editor `lb' .5 `lb' `top'
gr_edit .plotregion1.AddLine added_lines editor `lb' .5 `=`lb'+ 0.02' .5
gr_edit .plotregion1.AddLine added_lines editor `lb' `top' `=`lb'+ 0.02' `top'
gr_edit .plotregion1.AddLine added_lines editor `absent_itt_ub' .5 `absent_itt_ub' `top'
gr_edit .plotregion1.AddLine added_lines editor `absent_itt_ub' .5 `=`absent_itt_ub'- 0.02' .5
gr_edit .plotregion1.AddLine added_lines editor `absent_itt_ub' `top' `=`absent_itt_ub'- 0.02' `top'
gr_edit .plotregion1.added_lines[1].style.editstyle  linestyle( width(thin) color(red))
gr_edit .legend.plotregion1.key[4].DragBy 0 -26
gr_edit .legend.plotregion1.key[7].DragBy 0 -26
gr_edit .legend.plotregion1.label[12].DragBy 0 -26
gr_edit .legend.plotregion1.label[16].DragBy 0 -26
gr_edit .legend.plotregion1.label[9].DragBy 0 2.8
gr_edit .legend.plotregion1.label[13].DragBy .0 2.8
gr_edit .legend.plotregion1.label[6].DragBy 5.188841396576663 -4.425774722311327
gr_edit .legend.plotregion1.label[7].DragBy 5.188841396576663 -4.425774722311327
gr_edit .legend.plotregion1.label[6].style.editstyle color(`dkorange') editcopy
gr_edit .legend.plotregion1.label[7].style.editstyle color(red) editcopy
gr_edit .legend.plotregion1.label[6].style.editstyle size(small) editcopy
gr_edit .plotregion1.textbox1.style.editstyle vertical(top) editcopy
gr_edit .plotregion1.textbox1.style.editstyle size(small) editcopy

graph export "$Wellness_WhatDoesWWDo/results/figures/lit_review_absent_ITT.pdf", replace


twoway scatter counter bin if !out & order == 3, mcolor(`ltblue') msymbol(S) || ///
scatter counter bin if !out & order == 2, mcolor(`ltblue') msymbol(T) || /// 
scatter counter bin if !out & order == 1, mcolor(`ltblue') msymbol(O) || /// 
scatter counter bin if out & order == 3, mcolor(`ltblue') msymbol(Sh) || ///
scatter counter bin if out & order == 2, mcolor(`ltblue') msymbol(Th) || ///
scatter counter bin if out & order == 1, mcolor(`ltblue') msymbol(Oh) || ///
scatter counter bin if order == 4, msymbol(none) ||, ///
text(`itt_y' `itt_x' "★", color(`dkorange')) ///
yscale(off r(1 20)) ylabel(none,nogrid) graphregion(color(white)) legend(off) ///
xlabel(-.75 "-75%" -.5 "-50%" -.25 "-25%" 0 "0%" .25 "25%") xtitle("Percent Change in Absenteeism (ITT)") ytitle("") ///
name(absent_itt2, replace)

sum counter
local top = r(max) + 1
local lb = `absent_itt_lb'

gr_edit .plotregion1.AddLine added_lines editor `lb' .65 `lb' `top'
gr_edit .plotregion1.AddLine added_lines editor `lb' .65 `=`lb'+ 0.02' .65
gr_edit .plotregion1.AddLine added_lines editor `lb' `top' `=`lb'+ 0.02' `top'
gr_edit .plotregion1.AddLine added_lines editor `absent_itt_ub' .65 `absent_itt_ub' `top'
gr_edit .plotregion1.AddLine added_lines editor `absent_itt_ub' .65 `=`absent_itt_ub'- 0.02' .65
gr_edit .plotregion1.AddLine added_lines editor `absent_itt_ub' `top' `=`absent_itt_ub'- 0.02' `top'
gr_edit .plotregion1.added_lines[1].style.editstyle  linestyle( width(thin) color(red))
gr_edit .plotregion1.textbox1.style.editstyle vertical(top) editcopy
gr_edit .plotregion1.textbox1.style.editstyle size(small) editcopy

graph export "$Wellness_WhatDoesWWDo/results/figures/lit_review_absent_ITT_no_legend.pdf", replace

** IV effects for absenteeism

* Load estimates and calculate point estimate and CI

estimates use "`estimate_dir'/IV_ztreat_Lasso_sickleave_0816_0717.ster"

use "`lasso_dir'/lasso_data_sickleave_0816_0717.dta", clear

* JR edit
sum sickleave_0816_0717 [aweight= sickdays_eligible_0816_0717] if in_lasso_IV == 1 & hra_c_nomiss == 1

local absent_iv = _b[rhs]/(r(mean) - _b[rhs])
local absent_iv_lb = (_b[rhs] -  invnormal(0.975)*_se[rhs])/(r(mean) - _b[rhs])
local absent_iv_ub = (_b[rhs] +  invnormal(0.975)*_se[rhs])/(r(mean) - _b[rhs])

* Load OLS results

estimates use "`estimate_dir'/OLS_Lasso_sickleave_0816_0717.ster"

use "`lasso_dir'/lasso_data_sickleave_0816_0717.dta", clear
sum sickleave_0816_0717 [aweight= sickdays_eligible_0816_0717] if in_lasso_OLS == 1 & hra_c_nomiss == 0

local absent_ols = _b[rhs]/(r(mean) - _b[rhs])

* Load estimates from prior studies

use "$Wellness_WhatDoesWWDo/data/proc/lit_review/lit_review.dta", clear
drop if drop_dup
keep if estimate_type == "tot":estimate_type & outcome == "absenteeism":outcome

* Create an observation for our point estimate

set obs `=_N+2'
replace title = "Paper Estimate IV" in `=_N - 1'
replace title = "Paper Estimate OLS" in `=_N'
replace effect_percent = `absent_iv' in `=_N - 1'
replace effect_percent = `absent_ols' in `=_N'
gen main_iv = _n == _N - 1
gen main_ols = _n == _N
gen main = _n >= _N - 1

* Put all point estimates into one of 50 bins

sum effect_percent, detail

local max_x = max(ceil(100*r(max))/100,`absent_iv_ub')
local delta = (`max_x' + 1)/50
gen bin = .

forval x = 1/50 {

	replace bin = `max_x' - `x'*`delta' + .5*`delta' if effect_percent < `max_x' - `=`x'-1'*`delta' & effect_percent >= `max_x' - `x'*`delta'

}

egen bin_code = group(bin), label

di "`absent_iv_lb' `absent_iv' `absent_iv_ub'"

* Label point estimates that are out(side) of 95% CI
* Do it by bin, then move estimates near the boundary to the next bin if they are "out"

bysort bin_code (effect_percent): gen out = effect_percent[1] > `absent_iv_ub' | effect_percent[_N] < `absent_iv_lb'

replace bin_code = bin_code - 1 if !out & effect_percent < `absent_iv_lb'
replace bin_code = bin_code + 1 if !out & effect_percent > `absent_iv_ub'
replace out = 1 if out == 0 & (effect_percent > `absent_iv_ub' | effect_percent < `absent_iv_lb')

rename bin bin2
decode bin_code, gen(bin)
destring bin, replace

gen order = .
replace order = 5 if main_iv
replace order = 4 if main_ols
replace order = 3 if !main & study_type == "RCT":study_type
replace order = 2 if !main & study_type != "RCT":study_type & pre_post
replace order = 1 if !main & study_type != "RCT":study_type & !pre_post

bysort bin_code (order), sort: gen counter = _N + 1 - _n

sum bin if main_iv, meanonly
local iv_x = r(mean)
sum counter if main_iv, meanonly
local iv_y = r(mean)

sum bin if main_ols, meanonly
local ols_x = r(mean)
sum counter if main_ols, meanonly
local ols_y = r(mean)

twoway scatter counter bin if !out & order == 3, mcolor(`ltblue') msymbol(S) || ///
scatter counter bin if !out & order == 2, mcolor(`ltblue') msymbol(T) || /// 
scatter counter bin if !out & order == 1, mcolor(`ltblue') msymbol(O) || /// 
scatter counter bin if out & order == 3, mcolor(`ltblue') msymbol(Sh) || ///
scatter counter bin if out & order == 2, mcolor(`ltblue') msymbol(Th) || ///
scatter counter bin if out & order == 1, mcolor(`ltblue') msymbol(Oh) || ///
scatter counter bin if order == 4, msymbol(none) || ///
scatter counter bin if order == 5, msymbol(none) ||, ///
text(`iv_y' `iv_x' "★", color(`dkorange')) text(`ols_y' `ols_x' "☆", color(`dkorange')) ///
yscale(off r(1 20)) ylabel(none,nogrid) graphregion(color(white)) ////
xlabel(-.75 "-75%" -.5 "-50%" -.25 "-25%" 0 "0%" .25 "25%" .5 "50%" .75 "75%") xtitle("Percent Change in Absenteeism (TOT)") ytitle("") ////
legend(label(1 "RCT") label(2 "Pre/Post") label(3 "Other") ///
label(4 "RCT") label(5 "Pre/Post") label(6 "Other") ///
label(7 "OLS") label(8 "RCT") ///
order(- "IL Estimates:" 8 - "95% Confidence Interval" 7 /// 
- "Prior Studies:" - "★" - "[ ]" - "☆" ///
- "Not Ruled Out:" 1 2 3 /// 
- "Ruled Out:" 4 5 6) rows(4)) ///
name(absent_iv, replace)

sum counter
local top = r(max) + 1
local lb = `absent_iv_lb' - 0.02

gr_edit .plotregion1.AddLine added_lines editor `lb' .5 `lb' `top'
gr_edit .plotregion1.AddLine added_lines editor `lb' .5 `=`lb'+ 0.02' .5
gr_edit .plotregion1.AddLine added_lines editor `lb' `top' `=`lb'+ 0.02' `top'
gr_edit .plotregion1.AddLine added_lines editor `absent_iv_ub' .5 `absent_iv_ub' `top'
gr_edit .plotregion1.AddLine added_lines editor `absent_iv_ub' .5 `=`absent_iv_ub'- 0.02' .5
gr_edit .plotregion1.AddLine added_lines editor `absent_iv_ub' `top' `=`absent_iv_ub'- 0.02' `top'
gr_edit .plotregion1.added_lines[1].style.editstyle  linestyle( width(thin) color(red))
gr_edit .legend.plotregion1.label[6].DragBy 5.188841396576663 -4.425774722311327
gr_edit .legend.plotregion1.label[7].DragBy 5.188841396576663 -4.425774722311327
gr_edit .legend.plotregion1.label[8].DragBy 5.188841396576663 -4.425774722311327
gr_edit .legend.plotregion1.key[5].DragBy 0 -26
gr_edit .legend.plotregion1.key[8].DragBy 0 -26
gr_edit .legend.plotregion1.label[12].DragBy 0 -26
gr_edit .legend.plotregion1.label[16].DragBy 0 -26
gr_edit .legend.plotregion1.label[9].DragBy 0 2.8
gr_edit .legend.plotregion1.label[13].DragBy .0 2.8
gr_edit .legend.plotregion1.label[7].style.editstyle color(red) editcopy
gr_edit .legend.plotregion1.label[6].style.editstyle color(`dkorange') editcopy
gr_edit .legend.plotregion1.label[8].style.editstyle color(`dkorange') editcopy
gr_edit .legend.plotregion1.label[6].style.editstyle size(small) editcopy
gr_edit .legend.plotregion1.label[8].style.editstyle size(small) editcopy
gr_edit .plotregion1.textbox1.style.editstyle vertical(top) editcopy
gr_edit .plotregion1.textbox1.style.editstyle size(small) editcopy
gr_edit .plotregion1.textbox2.style.editstyle vertical(top) editcopy
gr_edit .plotregion1.textbox2.style.editstyle size(small) editcopy

graph export "$Wellness_WhatDoesWWDo/results/figures/lit_review_absent_TOT.pdf", replace

twoway scatter counter bin if !out & order == 3, mcolor(`ltblue') msymbol(S) || ///
scatter counter bin if !out & order == 2, mcolor(`ltblue') msymbol(T) || ///
scatter counter bin if !out & order == 1, mcolor(`ltblue') msymbol(O) || /// 
scatter counter bin if out & order == 3, mcolor(`ltblue') msymbol(Sh) || ///
scatter counter bin if out & order == 2, mcolor(`ltblue') msymbol(Th) || ///
scatter counter bin if out & order == 1, mcolor(`ltblue') msymbol(Oh) || ///
scatter counter bin if order == 4, msymbol(none) || ///
scatter counter bin if order == 5, msymbol(none) ||, ///
text(`iv_y' `iv_x' "★", color(`dkorange')) text(`ols_y' `ols_x' "☆", color(`dkorange')) ///
yscale(off r(1 20)) ylabel(none,nogrid) graphregion(color(white)) legend(off) ////
xlabel(-.75 "-75%" -.5 "-50%" -.25 "-25%" 0 "0%" .25 "25%" .5 "50%" .75 "75%") xtitle("Percent Change in Absenteeism (TOT)") ytitle("") ////
name(absent_iv2, replace)

sum counter
local top = r(max) + 1
local lb = `absent_iv_lb' - 0.02

gr_edit .plotregion1.AddLine added_lines editor `lb' .65 `lb' `top'
gr_edit .plotregion1.AddLine added_lines editor `lb' .65 `=`lb'+ 0.02' .65
gr_edit .plotregion1.AddLine added_lines editor `lb' `top' `=`lb'+ 0.02' `top'
gr_edit .plotregion1.AddLine added_lines editor `absent_iv_ub' .65 `absent_iv_ub' `top'
gr_edit .plotregion1.AddLine added_lines editor `absent_iv_ub' .65 `=`absent_iv_ub'- 0.02' .65
gr_edit .plotregion1.AddLine added_lines editor `absent_iv_ub' `top' `=`absent_iv_ub'- 0.02' `top'
gr_edit .plotregion1.added_lines[1].style.editstyle  linestyle( width(thin) color(red))
gr_edit .plotregion1.textbox1.style.editstyle vertical(top) editcopy
gr_edit .plotregion1.textbox1.style.editstyle size(small) editcopy
gr_edit .plotregion1.textbox2.style.editstyle vertical(top) editcopy
gr_edit .plotregion1.textbox2.style.editstyle size(small) editcopy

graph export "$Wellness_WhatDoesWWDo/results/figures/lit_review_absent_TOT_no_legend.pdf", replace

******************************
** Publication Bias Figures **
******************************

local dkorange = `" "232 74 39" "'
local ltorange = `" "244 165 147" "'
local dkblue = `" "19 41 75" "'
local ltblue = `" "137 148 165" "'
local dkgray = `" "gs12" "'
local ltgray = `" "gs12*.5" "'
local vltgray = `" "gs15" "'

graph set window fontface "Times New Roman"

use "$Wellness_WhatDoesWWDo/data/proc/lit_review/lit_review.dta", clear

local trim = 1000

sort pub_effect

drop if abs(pub_effect) > `trim'

** Sinificant studies

gen sig = abs(pub_effect)/pub_se >= 1.96

** generate data for z = 1.96 line

set obs `=_N + 3'

gen sig_y = .

replace sig_y = `trim'/1.96 in `=_N - 2'
replace sig_y = 0 in `=_N - 1'
replace sig_y = `trim'/1.96 in `=_N'

replace pub_effect = -`trim' in `=_N - 2'
replace pub_effect = 0 in `=_N - 1'
replace pub_effect = `trim' in `=_N'

gen sig_x = pub_effect

sort pub_effect

drop if pub_se > `trim'/1.96 & pub_se != .
keep if has_se
drop if drop_dup == 1

twoway line sig_y sig_x, lcolor(`dkgray') || /// 
	scatter pub_se pub_effect if !sig, msymbol(o) mcolor(`ltgray') || ///
	scatter pub_se pub_effect if sig, msymbol(o) mcolor(`dkblue') ||, ///
	graphregion(color(white)) ylabel(, angle(0) nogrid) xtitle("{&mu}") ///
	ytitle("{&sigma}") legend(label(2 "Not Significant") /// 
	label(3 "Significant") order(3 2)) name(pub_funnel, replace)
	
graph export "$Wellness_WhatDoesWWDo/results/figures/pub_bias_funnel.pdf", replace

** Make histogram of z scores

use "$Wellness_WhatDoesWWDo/data/proc/lit_review/lit_review.dta", clear

keep if has_se
drop if drop_dup == 1

gen new_z = pub_effect/pub_se

twoway histogram new_z, width(.2) density color(`dkblue') || ///
	pci 0 -1.96 .65 -1.96, lcolor(red) || ///
	, ylabel(0 "0" .2 "0.2" .4 "0.4" .6 "0.6", angle(0) glcolor(`vltgray')) ///
	graphregion(color(white) margin(r+5)) ytitle("Density") /// 
	xtitle("{&mu}/{&sigma}") text(.65 -4.5 "{&mu}/{&sigma} = -1.96", color(red)) legend(off) ///
	name(pub_z, replace)
		
graph export "$Wellness_WhatDoesWWDo/results/figures/pub_bias_z.pdf", replace

** EOF
